Shape of primordial non-Gaussianity and the CMB bispectrum 



J.R. Fergusson and E.P.S. Shellard 

Centre for Theoretical Cosmology, 
Department of Applied Mathematics and Theoretical Physics, 
University of Cambridge, Wilberforce Road, Cambridge CB3 OWA, United Kingdom 

(Dated: August 12, 2009) 

We present a set of formalisms for comparing, evolving and constraining primordial non-Gaussian 
models through the CMB bispectrum. We describe improved methods for efficient computation of 
the full CMB bispectrum for any general (non-separable) primordial bispectrum, incorporating a 
flat sky approximation and a new cubic interpolation. We review all the primordial non- Gaussian 
models in the present literature and calculate the CMB bispectrum up to / < 2000 for each different 
model. This allows us to determine the observational independence of these models by calculating 
the cross- correlation of their CMB bispectra. We are able to identify several distinct classes of 
primordial shapes - including equilateral, local, warm, flat and feature (non-scale invariant) - which 
should be distinguishable given a significant detection of CMB non-Gaussianity. We demonstrate 
that a simple shape correlator provides a fast and reliable method for determining whether or not 
CMB shapes are well correlated. We use an eigenmode decomposition of the primordial shape to 
characterise and understand model independence. Finally, we advocate a standardised normalisation 
method for /nl based on the shape autocorrelator, so that observational limits and errors A/tvl 
can be consistently compared for different models. 



INTRODUCTION 



Constraints on non-Gaussianity arguably provide the most stringent observational tests of the simplest inflationary 
paradigm and, in the near future, these limits are set to tighten substantially. Single-field slow-roll inflation predicts 
to high precision that the CMB will be a Gaussian random field and hence can be completely described by its angular 
power spectrum. However, if there was mechanism for generating some non-Gaussianity in the initial perturbations 
then its measurement would open up a wealth of extra information about the physics governing the early universe. 
Motivated by the first discussions of the local case which is dominated by squeezed states, non-Gaussianity is usually 
parametrised by /jvx, a quantity that can be extracted from CMB observations. The purpose of this paper is 
to apply and improve the methods detailed in ref. [1] for calculating the CMB bispectrum for any general (non- 
separable) primordial non-Gaussian model and, on this basis, to determine the extent to which competing models are 
independent and can be constrained by present limits on /nl- 

At present, CMB observations directly constrain only three separable primordial non- Gaussian models because 
of the calculational difficulties of using bispectrum estimators in the general case. These are the local model and 
the equilateral model (a separable approximation to the DBI inflation) and warm inflation with the latest (model- 
dependent) CMB constraints on f^L becoming [2-4] 

_ 4 <80, (1) 

-151</X< 253, (2) 
-375 < f% a L rm < 37 (3) 

The above are 2<j-limits on non-Gaussianity, though we note that there is a tentative claim [5] of a detection of /jy£ ai 
at almost 3<r, with limits 27 < /jy£ ai < 147 (2a). Given the growing range of theoretical possibilities for other (non- 
separable) non-Gaussian bispectra, an important goal is the development of methods which can be used to place direct 
constraints in the general case [6]. However, our purpose here is to note the models which are tightly constrained by 
the present limits and, conversely, to identify those for which further investigation is warranted. 

The bispectrum has been shown to be optimal for detecting primordial non-Gaussianity [7], but alternative ap- 
proaches seem to be able to produce comparable limits. A wavelet analysis of the WMAP 5 year data yields 
—8 < Jnl < HI (2a) [8] and positional information from wavelets can be used to examine the likelihood of specific 
features. Minkowski functionals provide a geometric characterisation of the temperature fluctuations in the CMB, 
yielding a slightly weaker constraint —70 < /nl < 91 (2a) [9]. Of course, the effects of non-Gaussianity are not only 
felt in the CMB but could be detectable in a wide range of astrophysical measurements, such as cluster abundances 
and the large scale clustering of highly biased tracers. However, given the imminent launch of the Planck satellite 
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with a projected constraint \/nl\ < 5 (2a) [10] (almost at the level of cosmic variance, \/nl\ < 3 (2<j)), the focus of 
this paper remains on the CMB bispectrum. 

In section 2, we review the relationship between the primordial bispectrum and its counterpart in the CMB, 
noting the importance of separability for present estimators. We also present a new analytic solution for the CMB 
bispectrum from a constant primordial model, against which it is useful to normalise other models. In section 3, 
we introduce the shape function relevant for nearly scale-invariant bispectra, and we review the current literature 
classifying current non- Gaussian models into centre-, corner- and edge- weighted categories, as well as those which 
are not scale-invariant. We use a shape correlator to forecast the cross-correlation of the CMB bispectra for all 
these different models, identifying five distinct classes of shapes. This work takes forward earlier discussions of the 
non-Gaussian shape (see, for example, [11-13]) but here we are able to directly compare to the actual CMB bispectra 
for all the models. Moreover, we propose a specific eigenfunction decomposition of the shape function which offers 
insight as to why particular shapes are related or otherwise. 

In section 4, we describe important improvements to the numerical methods which we use to make accurate 
computations of the CMB bispectrum for all the models surveyed [1]. The most important of these are the flat sky 
approximation and a cubic interpolation scheme for the tetrahedral domain of allowed multipoles. What was previously 
regarded as an insurmountable computational problem has now become tractable, irrespective of separability. In 
section 5, we present the main results detailing the cross-correlations for all the different models, while confirming the 
different classes of independent shapes previously identified, some of which remain to be fully constrained (even by 
present data). We compare the results of the shape and CMB bispectra correlators, noting the efficacy of the shape 
approach in identifying models for which quantitative CMB analysis is required. Finally, in section 6, we propose an 
alternative normalisation procedure for f^L which brings the constraints into a more consistent pattern, allowing for 
a model-independent comparison of the true level of non-Gaussianity. 



RELATING THE PRIMORDIAL AND CMB BISPECTRUM 



The primordial gravitational potential 3>(k) induces CMB temperature anisotropics which we represent using a^ m 's, 
that is, 

— = ^ai m Yi m {h) . 

Im 

The linear evolution which relates them is mediated by the transfer functions A/(fc) through the integral, 

— A,(fc)<£(k)Yi m (k). (4) 
The CMB bispectrum is the three point correlator of the a^ m , 

-Brnirn 2 m 3 = ( fl hmi fl l2m 2 %m 3 ) 5 (5) 

and so, substituting (4), we obtain 

($(k 1 )$(k 2 )$(k3))F iimi (ki)11 2m2 (k 2 )F i3 ™3(k3). (6) 

The primordial bispectrum is defined as 

<$(ki)#(k 2 )#(k3)> = (27r) 3 B $ (fc 1) k 2 , h) <y(ki + k 2 + k 3 ) , (7) 

where the delta function enforces the triangle condition. We replace the delta function with its integral form and 
expand the exponential into spherical harmonics. If we substitute this into equation (6) and integrate out the angular 
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parts of the three ki integrals, which yield delta functions, then we can remove the summation to obtain, 

£mim 2 m 3 = (j^j J dxdk^dk^xk^h) 2 'B*(k u k 2 , k 3 ) A^) A h (k 2 ) A h (k 3 ) 

ji 1 (k 1 x)ji 2 (k 2 x)ji 3 (k 3 x) J ^^ imi (x)^ 2m2 (x)^ 3m3 (x) . (8) 
The integral over the angular part of x is known as the Gaunt integral and has a geometric solution, 

Srnirn 2 rn 3 = J ^^l^i (x)lz 2 m 2 (x)^Z 3 m 3 (x) 



(2/i + 1)(2Z 2 + 1)(2Z 3 + 1) ( h h h \ ( h h I 



47r \ J \ mi rri2 rri3 



^3 



(9) 



I h h h \ 

where [ m m m J 1S ^he Wigner 3j symbol. First, as we expect the bispectrum to be isotropic, it is common to 
work with the angle averaged bispectrum, 

^iWa = ( m x m 2 m 3 ) ( a hmi^hm 2 a hm 3 ) • ( 10 ) 

We will find it more convenient to work with the reduced bispectrum where we drop the geometric factors associated 
with the Gaunt integral, 

ghhh = 0*1*2*3 ^ (11) 

In removing the 3j symbols it is important to remember some of the constraints that they place on the 1 values. The 
first constraint is that the sum of the three U must be even. The other is that if we treat the three U as lengths, they 
must be able to form a triangle. This is analogous to the constraint that the three ki are able to form a triangle, which 
is demanded by the delta function in (7), and in /-space is enforced through the x integral over the three spherical 
Bessel functions, which evaluates to zero if the corresponding triangle condition is violated. The reduced bispectrum 
is of the form, 

&Z1Z2Z3 = (j^j J dxdk 1 dk 2 dk 3 {xkik 2 k 3 ) 2 B^(k 1 ,k 2 ,k 3 ) 

A h (fex) A h (k 2 )A h (fc 3 ) j h (k lX )ji 2 (k 2 x)j h (k 3 x) . (12) 

In looking for a way in which we can simplify the integral associated with 6^/3, the most obvious place to start 
is to try to find simplifications for the primordial bispectrum B$. The most common, and simplest, approach is that 
pioneered by Komatsu and Spergel in [10]. Here they expand the primordial gravitational potential perturbation as 
a Taylor expansion around the Gaussian part, 

$(x) = $ G (x) + f NL (<t> G (x) - ($ G (x)) 2 ) , (13) 

where /nl parametrises the level of non-Gaussianity. Simple calculation results in a primordial bispectrum of the 
form, 

B*{k u fc 2 , k s ) = 2f NL (P*(fci)P*(fc 2 ) + P$(k 2 )P<s>(k 3 ) + P*(fc 3 )P*(fci)) • ( 14 ) 

This is known as the local model as the non-Gaussianity of the gravitational potential is local in space. Substituting 
this into equation (12) for the reduced bispectrum we see that, as the primordial bispectrum only consists of products 
of functions of a single ki, the integral can be separated into an integral over x of the products of integrals over /c, 

b hhh = 2f l £ c L al J x 2 dx (a h (x)(3 h (x)(3i 3 (x) + 2 permutations) , (15) 
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where, 



a i( x ) = \ f k 2 dk&i(k)ji(kx), 
h( x ) = \ J k 2 dkP^(k)Ai(k)ji(kx). 



(16) 
(17) 



This reduces the dimension of integration from four to two, and the separation allows us to calculate the reduced 
bispectrum easily. 

If we choose to work in the Sachs- Wolfe approximation, where we replace the transfer function with a Bessel 
function, 



(k) = ^ji ((r G - r dec ) k) , 



(18) 



the integral for the reduced bispectrum can be expressed in closed form and it is possible to derive an analytic solution 
in simple cases. The exact analytic result for the local model on large angles is 



where G(Zi , Z2, Z3) is shorthand for, 
G(7i, Z2, h) 



1 



1 



1 



h(h + i)Hh + 1) Hh + i)h(h + 1) h(h + i)Zi(Zi + 1) 



(19) 



(20) 



Here, we present a second analytic solution, that of the simplest primordial bispectrum that scales in the correct 
manner, 



B<s>(k 1 ,k 2 ,k 3 ) = l/(k 1 k 2 k 3 ) 2 , 
which we will denote as the constant model. Here the reduced bispectrum integral becomes, 



bhhh = Inl {j^j J x2dx n i=i4 (°> x ) ' 



Ii(p,x)= / k p dkji(k)ji(xk) . 



where 



The li(0,x) integral has a nice solution (see ref [14] p405), 



x > 1 



7T X 



X < 1 



7T X 



2 2/ + 1 ' ~ ^ ^ ' 2 21 + 1 

from which we can obtain an exact large-angle solution for the constant model bispectrum, 

h 



(21) 



(22) 



(23) 



(24) 



INL 



1 



27 (2Zi + l)(2Z 2 + l)(2Z 3 + l) 

3 



Pl POO 

J x h+l2+l3+2 dx + J x'^-^-^dx 



!. 



NL 



D{hMM); 



where we have defined, 

D{l u l 2 ,h) 



(2Zi + 1)(2Z 2 + 1)(2Z 3 + 1) 



li + h + h + 3 h + h + h 



(25) 



(26) 



Unfortunately the bispectrum signal is too weak for us to be able to measure individual multipoles directly from 
data, so to compare theory with observations we must use an estimator which sums over all multipoles. At the most 



5 



basic level estimators can be thought of as performing a least squares fit of the bispectrum predicted by theory, 
{ai irni ai 2rn2 ai 3Tn3 ) , to the bispectrum obtained from observations, ^ mi ^i 2 b m2 ^il m3 - If we ignore the effects of sky cuts 
and inhomogeneous noise the estimator can be written, 

c _ ^ { a hmi a l2m2 a hm 3 ) „obs n obs n obs / 97 \ 

Af 2 ^ C h C l2 C h lin " l3m3 ' { t) 

where 

M{B) 



B 2 

\ ^ hhh (28) 
\ ^ CiiCl 2 Ch 



The above estimator has been shown to be optimal [7] for general bispectra in the limit where the non-Gaussianity is 
small and the observed map is free of instrument noise and foreground contamination. This is of course an idealised 
case and we need to consider taking into account the effect of sky cuts, inhomogeneous noise, and beam effects which 
were considered in some detail in refs [12, 15]. 

Here, for clarity we consider only the simple form of the estimator (27), since all statements made about it can be 
extended to the more complete version. If we consider the local model we can rewrite it as, 

-1 n obs n obs n obs 

f = _LV nhhh h, 7 7 l im 1 a l 2 m 2 a l3m3 (9q \ 
J\f2 / j ^mim,2m3 U hhh C C C ' 

I i Tfi% 

From this we see that we only need to be able to calculate the reduced bispectrum, bi ± i 2 i 3 , from theory rather than 
the full bispectrum, {ai irni ai 2rn2 ai 3m , 3 ) , which would be much more challenging. For the local model we can use the 
separation of 6/ 1 / 2 / 3 , from equation (15), to separate the sums in the estimator, so, 

S = j^J d 3 xA( X ) (B(x)) 2 , (30) 

where, 

n obs 

^x^^Ta^-g^Wx) (31) 

Im 

n obs 

^(xlE^fifxl^fx). (32) 

Im 

From this we can conclude that if the primordial bispectrum is separable then we can overcome both the issue with the 
multi-dimensional integration, and the calculation of the 3j symbols. Thus separability has become the cornerstone 
of all non-Gaussian analysis. 



THE SHAPE OF PRIMORDIAL BISPECTRA 



Shape function 

Here we will introduce the shape function for the primordial bispectrum. As the power spectrum is constrained to 
be very close to scale invariant we expect that the bispectrum will show similar behaviour. Exact scale invariance for 
the local model results in an equal k primordial bispectrum of the form, 

Bi ocal (k,k,k)=6f NL ^-. (33) 

This equal-/c behaviour with B$(k, /c, k) oc k~ 6 turns out to be the expected scaling of a large number of non-Gaussian 
models, and so the difference between these models is only due to the dependence of the primordial bispectrum on 
the ratios k\ : k<i and k\ : ks [11]. As we already have the factor, (/C1&2&3) 2 , in the integral for the reduced bispectrum 
(12), it is natural to use it to divide out the scale dependence of the primordial bispectrum. Thus, we can define the 
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momentum dependence through a shape function S as, 

S(fci,fc 2 ,fc 3 ) = ^(k 1 k 2 k 3 ) 2 B^(k u k 2 ,k 3 ), (34) 

where TV is an appropriate normalisation, often taken to be TV = 1//nl- (Inconsistent definitions for /nl for different 
models mean that it is difficult to compare their observational limits; we will discuss this further in section 6, but 
here we will generally factor out N) . 

The two most commonly discussed models are the local model, 

S local (k u k 2 ,k 3 )oc^, (35) 



and the equilateral model 



^(fc 1; fc 2 ,fc 3 )oc^3. (36) 

iMll 



However, we should also keep in mind the constant model S const (ki,k 2 ,k 3 ) = 1, for which we have a large-angle 
analytic solution l 2l h)- Here in eqns (35-36), and throughout this section, we will adopt a shorthand notation 
for the possible combinations of wavenumbers that can contribute to the bispectrum (i.e. the simplest terms consistent 
with its symmetries): 

K p = ^2(ki) p with K = K X (37) 

i 

= 7^E(W fc i) 9 ( 38 ) 



K pq 



±pq 



Kpar = ^~ W (39) 

Pqr 

k ip = K P - 2(kif with h = k n • (40) 

where A pq = 1 + S pq and A pqr = A pq (A qr + 5 pr ) (no summation) . This notation significantly compresses the 
increasingly complex bispectrum expressions quoted in the literature. 

One of the first models discussed with a specific shape was Maldacena's calculation for single field slow-roll inflation 
[16]: 

S MaW (fei, fe, fc 3 ) cc (36 - 2 V )^- + e(K 12 + 8^) 

Am A 

« (6e - 2 V )S local (k 1 ,k 2 , k 3 ) + |eS e ^(fci, fc 2 , fc 3 )) , (41) 

where 5 ZocaZ and S e<?ni are normalised so that S local (k , fc , k) = S equi (k,k,k). While we know the predicted non- 
Gaussianity in this case is negligible, there are more recent models which yield similar combinations of equilateral 
and local terms which are measurable (e.g. non-local inflation [17]). We need to know, therefore, the extent to which 
we can distinguish between the relative contributions from these different shapes and the degree to which they are 
observationally independent. 



Shape correlator 

One obvious way to distinguish between models is to use the estimator discussed previously (27), replacing the 
observed bispectrum with one calculated from a competing theory, 

C(B ' B) - MbWW)^ c h c l2 c h • (42) 
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If the observational data contained a bispectrum of the form B[ il<2 ^ then C(B,B') is an estimate of the proportion of 
the correct f' NL we would recover by using an estimator based on Bi ± i 2 i 3 . However, this Fisher matrix approach is 
extremely computationally demanding as we must calculate the full bispectrum for each model before we can make 
any comparison. What we would like in addition, therefore, is a simple method allowing us to predict the value of the 
correlator directly from the shape functions, thus indicating cases in which a full Fisher matrix analysis is warranted. 

If we return to equation (12) for the reduced bispectrum and substitute the expression for the shape function we 
have, 

b h i 2 i 3 = Inl 0) jf <n> fe S(fci,fc2,fc 3 )A^ (43) 

where Vk is the area inside the cube [0, k ma x\ allowed by the triangle condition (refer to figure 1). The integral I G is 
given by, 

^Mtf-/.^«MMMMM- (44) 

So S(ki,k2,ks) is the signal that is evolved via the transfer functions to give the bispectrum today, with I G giving 
an additional, purely geometrical, factor. Essentially, I G acts like a window function on all the shapes as it projects 
from k to /-space, that is, it tends to smear out their sharper distinguishing features, but only erasing significant 
differences in extreme cases (as we shall discuss later). This means that the shape function S(ki, fe, £3), especially 
in the scale-invariant case, can be thought of as the primordial counterpart of the reduced bispectrum bi ± i 2 i 3 before 
projection. 

To construct a shape correlator that predicts the value of (42) correctly we then should consider something of the 
form 

F(S,S')= [ S(k 1 MM)S\k 1 MM)^ 1 MM)dV k , (45) 

where uo is an appropriate weight function. We note that authours in ref. [11] define a cosine between shape functions 
to be used as a correlator. However, we comment later on the quantitative differences of their definition, notably its 
breakdown in the non-scale invariant case. 

The question now is what weight function should we choose? Our goal is to choose S 2 uj in /c-space such that 
it produces the same scaling as the estimator B 2 /C 3 in /-space. Let us consider the simplest case where both 
k\ = &2 = ks = k and l\ = I2 = I3 = I. For primordial bispectra which are scale invariant, then, 

S 2 (k, k, k)u(k, k, k) cx u(k, k, k) . (46) 

If we work in the large angle approximation, and assume that I + 1 ~ Z, then we know C\ ex l/l 2 and from the analytic 
solutions G, D, (20, 26) that bin cx l/l 4 . The angle averaged bispectrum is related to the reduced bispectrum by, 



/ (2Z 1 + l)(2k + l)(2Z3 + l) h h h , , ) J7i 



For equal /, we can deduce that, 



B »i i )' ' 



The Wigner 3J symbol has an exact solution for which 



* (-IW^—JWWW « (-!)«/» JJ-I (49) 

oooy 1 > V3I+TV 3l\ (I/2)! 3 1 } VVsni' [ ' 
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0.55 
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1.00 



Table I: Shape correlations 



with the last expression using Stirling's approximation, /! « V2ttI(1 / ' ef . Combining these results gives, 



Hocr 1 , (50) 

and so we find that we should choose a weight function w(k, /c, fc) ex This is a very simple approximation which 
ignores the cross-sectional weighting inherent in (42). Our analysis of the behaviour in the (a, /^-slices with a constant 
primordial shape function (S(a,(3) = const) shows that B 2 /C 3 is flat in the interior and then grows to a finite value 
on the boundary. However this variation is confined to be very close to the boundary and we choose to neglect this 
effect. As a result we take the explicit flat k~ x weighting: 

w(h,k 2 , fc 3 ) = -j — —7—r -j , (51) 
fei + k 2 + k 3 

Note that this weighting does not incorporate damping due to photon diffusion at large /, edge effects or smoothing 
due to the projection from k to /-space. These could be included using phenomenological window functions, but 
our purpose here is simplicity. In any case, the choice of the weight function may significantly improve forecasting 
accuracy, but it does not impact important qualitative insights. 

With this choice of weight (51), the primordial shape correlator from (45) then takes the form 

Qfg g'\ _ S ) /j^) 

V ' } y/F(S,S)F(S',S') ' 

This correlator is very simple to calculate and allows us to quickly categorise new models and their degree of inde- 
pendence. 

Before surveying model shapes currently in the literature, we note the utility of this correlator (52) with the simple 
shape examples discussed above. The local and equilateral shapes (35, 36) yield a 46% correlation. Thus, supposing 
the universe to have local non-Gaussianity, a highly significant observation using a local estimator would be expected 
to produce a (less) significant result also for an equilateral estimator. Nevertheless, 46% is a relatively low correlation 
and the equilateral and local shapes can be regarded as distinguishable in principle. What of Maldacena's shape with 
e ~ r] (41) (as, for example, in m 2 (j) 2 inflation)? If this were indeed observable, it would yield the rather striking 
result that it is 99.7% correlated with the local shape (35) (see also [11]) and that it is only 53% correlated with 
equilateral (36). As we shall see, such strong correspondences between models with apparently different shapes are 
rather common. In this case, we need to go to a fine-tuned regime near r] « 2.84e in order for (41) to have an equal 
86% correlation to both local and equilateral shapes; it is clearly generically in the local class of models. Table I 
provides a summary of the correlations between all the shapes we discuss in the next sections. 

Finally we note that the definition in ref. [11] of a cosine correlator has three weaknesses. The first is that it is 
calculated only on a 2D slice, k\ = const, through the tetrahedron (in contrast to our 3D integration over the k — const 
slices in figure 1). This choice may be tolerable for comparing shapes which scale in exactly the same way as each 
slice will give the same correlation. However, if the two models being compared have differing scale dependence, like 
for comparisons with the feature model, then it will produce poor results. Even in the case when comparing shapes 
which scale in exactly the same way there is a second more subtle issue which makes the cosine correlator unsuitable. 
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A" 




Figure 1: (a) The region of /c-space allowed by the triangle inequality, i.e., for which the primordial bispectrum is valid. The 
red lines are ki = fo, &3 = 0; ki — /C3, ki = 0; = fci, = and the allowed region is in yellow, (b) This area can be 
parametrised into slices represented by the green triangle and the distance 2|/c|/v / 3 of the centre of the triangle from the origin. 



The region covered in the bispectrum correlator is the intersection of the cube defined by [0, 1 max] and the tetrahedron 
denned by the triangle condition on the three U and so we should cover a similar region in /c-space. If we just look at 
the correlation on one slice then we miss the effect the shape of the region has on the result. If we think of the region 
as being composed of many parallel slices then some will be incomplete due to the effect of restricting the individual 
h < k max . Different slices will give different correlations depending on how much they have been cut and so no slice 
is truly representative of the true correlation. The third and related problem if that the weight in (51) is required to 
give accurate representation of the CMB correlation using shape functions in /c-space. 



Shape decomposition 

Given strong observational limits on the scalar tilt we expect all shape functions to exhibit behaviour close to 
scale- invar iance, so that S(ki,k2,ks) will only depend weakly on |fc|. Here, we choose to parametrise the magnitude 
of the fci's with both |fc| = (k\ + k\ + &1) 1//2 and the semi-perimeter, 

k = ^(h + h + ks). (53) 

A consequence of this scaling behaviour is that the form of the shape function on a cross-section is essentially 
independent of fc, so that for the models under consideration we can write 

S(k u k 2 ,ks)=f(k)S(kuk 2 ,k 3 ). (54) 

where 

7 kl 7 ^ 7 ^ /rr\ 

1 = T' 2 = T' 3 = T' (55) 

and we note that fci + k 2 + ks = 2. Since we are restricted to the region where the three ki are able to form a 
triangle by momentum conservation, we will reparametrise the allowed region to separate out the overall scale k from 
the behaviour on a cross-sectional slice Sk- This two-dimensional slice is spanned by the remaining coordinates (see 
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Figure 2: Coordinate transformation from a shaded subtriangle on the equilateral (a, (3) slice to the uniform square (x, y) 
domain suitable for an eigenmode expansion, i.e. a = 1 — x, (3 = yx/3. Here, we illustrate the transformation with contour 
plots for the warm inflation shape function. Note the non-trivial behaviour in the corner region where the function diverges 
and the sign changes. 

figure 1), 

fci = k (1 - f3) , 

k 3 = h(l-a + (3) . (56) 

The surface k = const, defines a plane with normal (1, 1, 1) at a distance 2|fc|/\/3 from the origin. Our new parameters 
a, /?, parametrise the position on the triangular domain formed by the intersection of the tetrahedral region and that 
plane [18]. They have the following domains, < /c < oo, < /? < 1 and — (1 — (3) < a < 1 — /3. In this parametrisation 
we can re-write shape function (54) and the volume element respectively as 

S(ki, &2, ks) = f(k)S(a, (3) , dVk = dk\dk2dks = k 2 dkdad(3 . (57) 

Here, we note that f(k) « const, for all the model shapes to be discussed in the next section, with the exception of 
the feature and oscillatory models. 

Having restricted our discussion to two-dimensional (a,/?)-slices, we now note that bispectrum symmetries are 
such that we need only characterise the shape on one sixth of this domain (refer to figure 2). This is a right-angled 
subtriangle with corners defined by the centre of the original triangle plus any corner together with the midpoint of 
an adjacent side. Here, we choose the bottom right triangle (containing fc 3 — > 0), with corners (~, 0), (1, 0) and (0, 
0) (i.e. the shaded region in figure 2). In order to set up a straightforward eigenfunction decomposition, we make the 
following coordinate transformation to take our subtriangle to a unit square, 

a = 1 — x, (3 = yx/3, (58) 

with < x < 1 and < y < 1. Analogous to polar coordinates at r = 0, this transformation expands the ks = 
corner (see figure 2) and our Sk volume element becomes 

da dp = x dx dy . (59) 

Here, we note that x and ks share the same squeezed limit (x — > as ks — > 0). 
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With the simple weight function w(x,y) = x we can now decompose an arbitrary shape function S(x,y) defined on 
Sk into a sum, 



S(x,y) = ^2c mn X m (x)Y n (y) , 



(60) 



consisting of products of orthogonal eigenfunctions X m and Y n defined on the unit interval. One possible choice would 
be Bessel functions X n (x) = J p (X pn x) (given the weight w = x) and trigonometric functions Y m (y) = Asm{ny) + 
Bcos(ny). However, the shape function in general is neither periodic nor vanishing on the boundary, leading to an 
ill-conditioned problem with poor series convergence. Given these arbitrary boundary conditions, a better choice 
employs Legendre polynomials P n (y) in the ^/-direction and analogous radial polynomials R m (x) in the x-direction. 
The domain < y < 1 requires shifted Legendre polynomials P n (y) which, unit normalised, become 



P (y) = l, Pi(ff)=V3(-l + 2y) P 2 (y) = V5(l-6y + 6y 2 ) 
The eigenfunctions in the x-direction can be created using the generating function, 



(61) 



R n (x) 



N 



1/2 


1/3 


1/4 . 


. l/(n + 2) 


1/3 


1/4 


1/5 . 


. l/(n + 3) 


1/n 


1/n+l 


1/n + 2 . 
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X 
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(62) 



with the first unit normalized polynomials given by, 

Ro(x) = y/2, Ri(x) = -4 + 6x . 



R 2 (x) = ^/§(?> - I2x + lOx 2 ) , 



(63) 



We can now find an eigenmode decomposition (60) on the domain Sk for any given shape function S(x, y) with the 
expansion coefficients given by 



ff 

Jo Jo 



Rm (x) P n (y) S(x,y)xdxdy . 



(64) 



Exploiting eigenmode orthogonality, the counterpart of the correlator (52) between two shapes 5, S f on the Sk slice 
then is 



C Sh (S,S')= f S(a,f3)S f (a,f3)dadf3 = J2 ( 

J S k n.m 



(65) 



where we have assumed unit normalised S, S' on Sk- Even for scale-invariant shapes, the slice correlator (65) is not 
identical to the overall shape correlator (52) because the integration domain for the latter includes the cubic region 
[0,k mas]) that is, for large k > k max it integrates over interior regions of the slices which weights the centre more 
heavily. (The integration domain is illustrated explicitly in figure 12.) Nevertheless, Cs k (S,S f ) is in close agreement 
with C(S, S') for highly correlated shapes and is able to reliably distinguish between independent shapes. 

Eigenmode expansion matrices (c mn ) illustrating the two qualitatively different results for the equilateral S equ u and 
local Siocal shapes are, respectively, 



(66) 



The eigenmode coefficients converge rapidly for the equilateral shape (left) and it can be very well approximated by 
just three linear terms (with which it is 98% correlated): 
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Sequii(x,y) « 0.88 + 0.35J?i(a:)+0.27Pi(y) = -0.74+ 1.65a: + 1.76y. 



(67) 
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Figure 3: Eigenvalues for the two-dimensional eigenmode expansion of the warm inflation shape function. Here, we denote the 
conventions with m (for Rm(x) incrementing in the horizontal x-direction and the n (for P n (y)) in the vertical ^/-direction. 
Note the dominance of R m (x)Pi(y) modes. The colour coding (used also in figure 4) is such that only blue and red colours can 
contribute at above the 10% level to the autocorrelator Ck(S,S f ), with yellow and pale green below 1%. 



The local shape, on the other hand, which is divergent except for a cut-off at k m i n /k max ~ 2/Z ma£C , oscillates as 
~ ±l/y/m in the R m (x) modes (for n = 0), although it converges rapidly for the higher P n (y) (n > 1). It is 
immediately apparent that the 46% correlation between local and equilateral shapes discussed earlier arises primarily 
from the dominant constant term cooCq in the expansion (65). Removing the constant mean from the local shape, 
then yields a small -15% anticorrelation between the models, thus suggesting possible strategies for distinguishing 
them. 

In the discussion that follows we shall illustrate the eigenvalues graphically as shown in figure 3 for the warm 
inflation shape. Like the local model, the warm shape is corner- weighted, however, the dominant presence of the 
R m (x)Pi(y) modes (see also figure 2), which are orthogonal to the R m (%) modes in the local case, implies that the 
two shapes exhibit little correlation (only 30%) and can be regarded as independent. Figure 4 provides a summary of 
the largest eigenvalues for all the shapes we discuss in the next sections. 



Equilateral triangles — centre-weighted models 

We begin this brief survey with the shape functions which are most easily characterised, those bispectra dominated 
by contributions from nearly equilateral triangle configurations, k\ « ki ~ ks. While these might be well-behaved 
shapes, they are not necessarily the best-motivated physically. Equilateral non-Gaussianity requires the amplification 
of nonlinear effects around the time modes exit the horizon, which is not possible in a slow-roll context for vanilla 
single field inflation. Instead, the kinetic terms in the effective action must be modified as in the Dirac-Born-Infeld 
(DBI) model [19] or by explicitly adding higher derivative terms, such as in K-inflation (see, for example, ref. [20]). 
The resulting corrections modify the sound speed c s , acting to slow the scalar field motion and, when the field theory 
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Figure 4: Eigenmodes for different shape functions using the conventions and colour scale defined in figure 3. Strong similarities 
are apparent between the equilateral family of models which are all highly correlated and would prove very difficult to distinguish 
observationally. The independence of the local and warm models is also apparent from the orthogonality of the dominant 
eigenmodes. 



is coupled to gravity, inflation is able to take place in steep potentials. For DBI inflation, this leads to non-Gaussianity 
being produced with a shape function of the form [19, 21] 

S DBI [k u fc 2 , h) (x - * (K& + 2K 14 - 3K 23 + 2Kns ~ 8^122) , (68) 
KmK 2 

where we have used the compact notation of equation (37). Another example of a model with non-standard kinetic 
terms is ghost inflation [22]. Here, a derivatively-coupled ghost scalar field <p is responsible for driving inflation. When 
(j) < the potential can be thought of as flat, but <\> 7^ and so the field continuously evolves towards = 0, where 
inflation ends. In this model the dominant effect for the perturbations comes from the trilinear term in the Lagrangian 
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Figure 5: The shape function of models in the equilateral class. Clockwise from top left we have the equilateral, DBI, single 
and ghost models. All four of these models have the majority of their signal concentrated in the equilateral limit corresponding 
to the centre of the triangle. Despite significant variations in the flattened limit, particularly around the edges of the triangle, 
all are strongly correlated by 96% or greater to the equilateral model 



which naturally leads to a non-zero bispectrum. The shape function for this model is of the form, 



S 9host (k 1 ,k 2 ,k 3 )oc2 



plus two permutations, where, 



Km 



-Re 



7] k 3 k 



32 



(69) 



(70) 



General non- Gaussian shapes arising from modifications to single field inflation have been extensively reviewed in 
ref. [20]. Using a Lagrangian that was an arbitrary function of the field and its first derivative, they were able to 
identify six distinct shapes describing the possible non- Gaussian contributions. Half of these had negligible amplitude 
being of the order of slow-roll parameters (two already given in (41)). Of the remaining three shapes [20] (see also 
[23]), one was believed to be subdominant, the second recovered the DBI shape (68), leaving a third distinct single 
field shape of the form, 



S single (h,k 2 ,k 3 ) 



oc 



Kill 

K 3 



(71) 



The sub-dominant term is a complex combination of special functions (somewhat like ghost inflation (69)) with inde- 
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Figure 6: The shape function of models in the local class. On the left is the usual local model, while on the right is the non-local 
model, which is virtually identical. These two models are very highly correlated and so fall into the same class 

terminate parameter values; we will neglect it in the subsequent discussion. Finally, we recall the original equilateral 
shape (36), noting that it was introduced not because of a fundamental physical motivation, but as a separable 
approximation to the DBI shape (68) [11]. 

Using the shape correlator (52) and the shape decomposition (64) introduced above, we can make a preliminary 
comparison of the four equilateral shapes — DBI, ghost, single and equilateral — which are illustrated on a k — const, 
slice in figure 5. The results of the shape correlators are given in the table I. Despite the apparent visual differences 
between these shapes, particularly near the edges of the triangular domain, there is at least a 96% or greater correlation 
of each to the equilateral shape (36). These particular centre- weighted shapes must then be regarded as a single class 
which would-be extremely difficult to distinguish observationally (see later for the CMB bispectrum discussion). 

The underlying similarity of the shapes is evident from the magnitudes of the eigenvalues illustrated in figure 4. 
Each shape has a dominant constant term coo and can be very well- approximated by a linear polynomial S(x, y) ~ 
coo + c\qR\(x) + Co i Pi (y) on the relevant subdomain (within a 96% correlation as previously (67) for the equilateral 
shape). The greatest difference occurs between the single and ghost shapes, with the former completely dominated 
by the coo and the latter having significant coi and en eigenvalues. Independent equilateral shapes are possible, in 
principle, but they must go further than the ghost shape in suppressing the constant term coo in favour of eigenmodes 
with greater internal structure (n, m > 1). 

Squeezed triangles — corner-weighted models 

The local shape covers a wide range of models where the non-Gaussianity is produced by local interactions. These 
models have their peak signal in "squeezed" states where one hi is much smaller than the other two, this is because 
non-Gaussianity is typically produced on superhorizon scales. The simplest case is that of single-field slow-roll inflation 
(41) [24], which as we have seen is dominated by the local shape. The non-linearities produced are tiny and /jy£ aZ is 
constrained to be of order slow-roll parameters [16, 24-26]. The production of non-Gaussianity during multiple field 
inflation [27-32] shows much greater promise (see, for example, recent work in refs. [33, 34] and references therein). 
Here non-Gaussianity is created by the inflaton when it follows a curved trajectory in phase space, during which 
isocurvature perturbations are converted into adiabatic perturbations [35, 36]. The magnitude of the non-Gaussianity 
generated is normally around f l § c ^ 1 ~ 0(1), which is at the limit for Planck detection. Significant f l § c ^ 1 can be 
produced in curvaton models [37-39] where the adiabatic density perturbation is generated after inflation by an 
initially isocurvature perturbation in a light scalar field, different from the inflaton. The non-Gaussianity generated 
in this scenario can be as large as, /jy£ aZ « 0(100). Large /jy£ aZ can be generated at the end of inflation from massless 
preheating or other reheating mechanisms [40-44]. After slow-roll inflation ends, the inflaton oscillates about its 
minimum and decays. Preheating occurs when a light field oscillates in resonance with it, taking energy from the 
inflaton, so its amplitude grows. The amplitude of the resonant field eventually becomes so large that its dynamics 
become non-linear and this non-linearity is transferred to the density perturbations. It is claimed this process can 
generate enormous non-Gaussianity, /jy£ aZ ~ 0(1000), which is already tightly constrained by observation. 
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Figure 7: The shape function of warm inflation before (left) and after smoothing (right). Note the pathological behaviour of 
the squeezed states in the corners for the unsmoothed model. The derivation for this shape breaks down when || < ^p^- as 
kl — > 0, where T is the friction coefficent due to dissipative effects. The smoothed model presents a more reasonable profile, 
but results are cut-off dependent. 



The local shape is strongly motivated because it appears in models that use standard kinetic terms in the action, 
smooth potentials without exotic couplings and which assume the standard Bunch-Davies vacuum. We note, however, 
that it also occurs in other contexts. Significant local non-Gaussianity can appear in models based on non-local field 
theory, such as p-adic inflation [17]. In these models slow roll inflation is again able to occur in very steep potentials. 
Like single field slow-roll inflation, the predicted non-local shape function is a combination of local (35) and equilateral- 
like (36) shape functions (see also refs [45-47] for its origin). However, the combination is even more heavily weighted 
than (41) towards the local shape (with the relative ratio given roughly by the number of e- foldings). Consequently, 
the non-local shape is (paradoxically) completely indistinguishable from the local shape (and is subsumed in this 
class henceforth). A comparison of the two shapes can be seen in figure (6). The ekpyrotic model can also generate 
significant /jy£ a ' [48-52]. Here the density perturbations are generated by a scalar field rolling in a negative exponential 
potential, so non-linear interactions are important, and large local non-Gaussianity can be produced, f l £f l « 0(100). 

Finally, we note that warm inflation scenarios, i.e. models in which dissipative effects play a dynamical role, 
are also predicted to produce significant non-Gaussianity [53-55]. Contributions are again dominated by squeezed 
configurations but with a different more complex shape, 

S warm (k u fc 2 , k s ) (x ^—(K 45 - K 27 + 2K 225 ) . (72) 

^333 

As we can see from figure 2, the squeezed limit contains an orthogonal sign change as the squeezed limit is approached 
h - 0. 

It is immediately apparent that we need to introduce a cut-off in order to normalize the squeezed shape functions 
for the correlator (52). This logarithmic divergence does not afflict the CMB bispectrum Bi ± i 2 i 3 because we do not 
consider contributions below the quadrupole 1 = 2. Given the cut-off at large wavenumbers where k ma x is related 
through a flat sky approximation to the largest multipole Z macc , we can similarly define k m i n w (2/l max ) k ma x- There is 
only a weak dependence on the precise value of k m i n . We note that a more serious concern is the apparent breakdown 
of the approximations used to calculate the warm inflation shape near the corners. In the absence of a specific 
prescription for this asymptotic regime, we have to explore the dependence of an edge cut-off and/or smoothing. We 
remove the divergence in the squeezed limit, k\ — > 0, by truncating the shape function when k\/(k 2 + ks) < 0.015. 
We smooth the resulting discontinuity by applying a Gaussian window function on the cross sectional slices with a 
FWHM of 0.03/(fci + k 2 + ks). The result of this applying this process to the shape function can be seen in figure 
(7), with the model denoted warmS. 

We note that the local and warm shape functions are only correlated at the 33% level, despite the primary contri- 
bution coming from squeezed states in both cases. As discussed previously, this is because their dominant eigenvalues 
c m o (local) and c m i (warm) correspond to orthogonal eigenmodes. 

The local shape is modestly correlated at the 40-55% level with the equilateral shapes (qualitatively in agreement 
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Figure 8: The shape function of the flattened model before (left) and after smoothing (right). The approximation diverges in 
the flattened limit, whereas it should oscillate and decay to zero. Here, we set the boundary to zero then apply a Gaussian 
smoothing to obtain a more reasonable profile. 



with [11], though not quantitatively presumably because of the different weighting). In contrast, however, the local 
contribution from the constant term coo ~ 0.55 is relatively small. Thus removing the coo term from the local 
estimator eliminates most of the equilateral correlations while leaving 70% of the local signal (i.e. in the autocorrelator 
Cs k (S, S)). Thus, we propose subtraction of the constant term as a significant test of the local model. 



Flattened triangles — edge-weighted models 

It is possible to consider inflationary vacuum states which are more general than the Bunch-Davies vacuum, such as 
an excited Gaussian (and Hadamard) state [56]. Observations of non-Gaussianity in this case might provide insight 
into trans-Planckian physics. The bispectrum contribution from early times is strongest for flattened triangles with, 
for example, ks « k\ + k<i- In the small sound speed limit c s « 1, the primordial bispectrum could be significant with 
a shape given by [20] 

Sf lat (k u k 2 ,k 3 ) cc -L_ (K 12 - K 3 ) + 4^§^ . (73) 

rC-^ K"2 

Unfortunately, as this analytic approximation diverges on the entire boundary of the allowed region, any integrals 
over the bispectrum are unbounded. In principle, however, this divergence should be cut-off near the edges, though 
the nature of the asymptotic behaviour which replaces it is poorly understood. In order to obtain even a qualitative 
picture for the flat shape function we must truncate it in some way. We follow the same procedure as for the warm 
model, removing the section ki/K < 0.03, then applying the same Gaussian filter to remove the discontinuity. We 
refer to this shape as S^ latS . Plots of the flattened model before and after smoothing can be seen in figure (8). 

Reflecting its distinctive properties, the flat shape is poorly correlated with most of the other shapes, with a 
particularly striking absence of any correlation with the orthogonal warm shape. Having a dual divergence (xy) -1 
means that the eigenvalues are spread more thinly and widely than the corner- weighted models with a smaller constant 
term. Nevertheless, the fiat shape would be most susceptible to confusion with the local shape with which it has a 
62% correlation. 



Features — scale-dependent models 

Finally, there are models in which the inflation potential has a feature, providing a break from scale-invariance 
and introducing large scale power where it is deemed to be indicated by observation. This can take the form of a 
either a step [57] or a small oscillation superimposed onto the potential [58]. Analytic forms for both these three 
point functions have been produced in [59]. However, these approximations are both somewhat simplistic and so 
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are unsuitable for a detailed analysis, other than as preliminary check of their correlation to other models. The two 
approximations are of the form, 



K 



Sf eat (k 1 ,k 2 ,k 3 )<xsin(- + P), 
S osci (k 1 ,k 2 ,h) oc sin(Cln(fO + P) , 



(74) 
(75) 



where k* is the associated scale of the feature in question, C is a constant and P is a phase factor. The correspondence 
of these analytic approximations to the full shape function can be seen in figure (9) 

Results for the shape correlator for a particular feature model (with k* « £*/ r o ano ^ ^* = 50), are given in table (I). 
It can be seen to be essentially independent of all the other shapes. Obviously this is because all variations in feature 
model occur in the if -direction which is orthogonal to the (a, /?)-slice - it only shares the constant term in common. 




Figure 9: The scaling of the shape functions of models with features in the potentials. These are plots of the scaling of the 
central value S(k, k, k) for a model with a single feature on the left, and a potential with an oscillatory component on the right. 
The blue dashed line is the correct numerical result, the red solid line is the simple approximation quoted earlier (these plots 
approximate those in [59]). 

We conclude, from this brief survey of the literature, that we can identify the feature model in a fifth distinct 
category beyond the equilateral, local, warm and flat shapes. We shall now turn to the much more formidable task 
of calculating the CMB correlators directly in order to determine the accuracy of our shape correlator analysis. 

CMB BISPECTRUM CALCULATION METHODOLOGY 
Numerical approach 

It is not feasible to directly evaluate the bispectrum for a completely general model. However, provided the shape 
function obeys a mild separability ansatz then the reduced bispectrum integral can be re-written in a tractable form. 
The method is based on the splitting of the shape function (34) into scale and scale-free parts (57), S(ki,k2,ks) = 
f(k)S(a, (3) , as discussed in the previous section, that is, an ansatz which applies to all the models under discussion. 
By using this decomposition with the reparameterisation into rescaled wavenumbers fci, fe, ks from (55), we can 
rewrite the integral for the reduced bispectrum (43) in a simple form 



Inl ( - 

7T 



^^^^(^^(^^^^ajA^^OA^^)^^) 
x l£i 2 i 3 (kk u kk 2 ,kk 3 ) 



^5(a,^)/^ 2 , 3 (a,/3)^ 2 , 3 (a,/3), 



(76) 
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where Sk is the cross-section spanned by a and (3 from (56) and dSk = dad/3. Here we have made the definitions, 

. . . f°° f(x] 
I? lhl3 (k 1 ,k 2 ,k 3 )= / dx J -^A ll (xk 1 )Ai 2 (xk 2 )Ai 3 (xk s ), (77) 
Jo x 

^2/3(^1^2,^3) = / dxx 2 j h (xk 1 )ji 2 (xk 2 )ji 3 (xk 3 ) . (78) 
Jo 

We refer to eqn (77) as the transfer integral and eqn (78) as the geometric integral. With this decomposition we 
have reduced the number of dimensions in the integral from four to three and we have also trapped all the highly 
oscillatory behaviour into the two one-dimensional integrals, I T and I G . Having achieved this, the remaining two- 
dimensional integral over Sk has very mild oscillatory behaviour and only requires a similar number of points as the 
one-dimensional integrals to evaluate accurately. 




Figure 10: A plot of the three spherical Bessel functions and the resulting integrand 1° . The top three plots are of 31 (ax) for 
different values of / and a, with (20,0.05) (top), (40,1.0) (second) and (50,0.95) (third). In the bottom diagram, the resulting 
integrand for I G is plotted. From this we can see that the peak sections of the two middle plots are ignored and it is their 
tails that are picked out in the integrand. Thus, for accurate evaluation we must calculate the transfer functions and Bessel 
functions over a much larger range than is required to evaluate the power spectrum. 

We will now detail the numerical methods used to evaluate both the one-dimensional transfer and geometric integrals 
and the two-dimensional integral over the triangular domain Sk- We use a modified form of the CMB bispectrum 
code already presented in [1], so here we will focus on substantial recent improvements. 

To evaluate the bispectrum using this formalism we must first be able to compute the two ID integrals, (77) and 
(78), for every combination of the three k{ possible in the triangular domain Sk- The transfer integral converges 
quickly, as l//c 4 , so while it is highly oscillatory it does not pose an enormous challenge. Also the transfer functions 
truncate at large values of k due to photon diffusion and so the integral naturally terminates. The geometric integral 
only converges slowly, as 1/fc, and so constitutes the majority of time in calculating the integrand for each point in 
the triangular domain Sk- 

To evaluate these two one-dimensional integrals we first need to obtain the transfer functions and Bessel functions 
that make up their integrands. We cannot simply output them from the currently available CMB temperature 
anisotropy codes, like CAMB and CMBFast, as their ranges are insufficient for our purposes. This is due to the 
rescaling of the functions by the three fc^, which range between [0,1]. Although, due to the constraint k\ + k 2 + k% = 2, 
only one can be small at a time. This has the effect of stretching the functions relative to each other. Both the transfer 
functions and Bessel functions have a similar form. They begin with a long section which is approximately zero before 
beginning oscillations which decay as k~ x . This means that they can pick out sections of the other functions that 
would have been unimportant for the calculation of the power spectrum. This effect can be clearly seen when we plot, 
in figure (10), the three individual functions ji(xk) and the integrand of the geometric integral I G for a point in the 
triangular domain Sk generated for the reduced bispectrum point 6204050- It is clear that it is the tail of the second 
and third Bessel functions which is important for calculation, rather than the initial region where their individual 
signals are largest. Thus, to accurately calculate the two ID integrals we need both the transfer functions, and Bessel 
functions, to cover a much larger range of fc, with a much better resolution, than is required to evaluate the power 
spectrum. 
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Our approach is to output the source function as calculated by one of the CMB anisotropy codes, having first 
increased the range and resolution of k and the resolution of r. We then calculate the transfer functions ourselves. 
The Bessel functions are calculated using a similar method to that in CMBFast. Both are then stored in tables for 
later use. This fixes us to a particular cosmological model unless we wish to keep regenerating the tables. As as 
we are primarily concerned with selecting primordial models, and the bispectrum is too small to usefully constrain 
cosmological parameters, this is a minor concern at this stage. 

At the beginning of a calculation, we read in the tables selecting the rows that correspond to the relevant Vs and 
interpolate them into a cubic spline. For each point in the triangular domain Sk we then take the three k{ and 
calculate the three function values corresponding to the rescaled points, kix. The three rescaled functions can then be 
multiplied together to form the integrand and we then use a cubic spline to evaluate the integral. The use of splines 
significantly reduces the resolution needed for accurate evaluation as compared to using simple linear interpolation 
(used previously [1]), usually by an order of magnitude. Developing faster integration methods, like the search for an 
analytic solution to the geometric integral, remains an interesting avenue for future work with scope for significant 
efficiency gains. With this method there are only two main parameters that control convergence, the resolution and 
the range. These have undergone extensive experimentation and minimal values for sub 1% accuracy have been found. 
The calculation for the 2D integral was completed using the same adaptive method used in [1]. 

This approach allows us to accurately calculate the bispectrum for a broad range of primordial non-Gaussian 
models. If we wish to determine the bispectrum at the resolution of Planck, l m ax 

« 2000, then the possible allowed I 

configurations require over 600 million integrations. Fortunately, there are several techniques we can use to make this 
problem tractable. These calculations naturally coarse-grain the computational work either through the sampling of 
ID integrals on the 2D triangular grid or, at a higher level, simply by evaluating the bispectrum at different multipole 
values. Secondly, the problem is well-suited to parallelisation on a large supercomputer or cluster and this has been 
achieved with the present code using a message passing interface (MPI) implementation which significantly reduces 
calculation time-scales. However, there are two further methods we use to dramatically speed up calculation which 
are detailed in the following subsections. 



Flat sky approximation 

If we calculate the reduced bispectrum when all three U are large then we are considering very small angles in the 
sky. If the angles are small then the curvature of the surface of last scattering is small and the so the sky can be 
approximated as flat. This allows us to greatly simplify the integral for the reduced bispectrum. We use the flat sky 
methodology which first appeared in [11] beginning by expanding the temperature perturbation into plane waves, 

e(x,n) = |^a(l)e- il - n a(l) = J <PnQ{x, n)e fl - n . (79) 

and so, 

o(l) = J £° dTS(k,T)e- ikZ< - T °- T) 5 2 (k"(ro - r) - l) , (80) 

where we have split k into the part parallel to the tangent plane, k", and the part perpendicular, k z . We now form 
the three point correlator for the flat sky a(l)'s, 

(a(li)a(l 2 )a(l 3 )) = J dr^d^d 3 k x d 3 k 2 d 3 k z 5 (J2 k i) S 2 (J2 k h B ^( k u k 2,h) 
5(fci,T 1 )5(/c2,T 2 )S , (fc3,T3)e- ifc i( T0 - Tl )e- ifc ^ T0 - T2 )e- ifc 3( T °- T3 ) 
S 2 (kf (to - n) - h) S 2 (k!(r - r 2 ) - 1 2 ) S 2 (kjj(r - r 3 ) - 1 3 ) . (81) 

To integrate out the three delta functions for k" we must assume that the variation in B$(ki, &2, /C3) is small in the 
k" direction across the width of last scattering. This allows us to use an average value for k" of k" = 1/(tq — tr). 
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After substitution this gives, 

<o(li)o(l 2 )o(l3)> = (to - t r ) 2 S 2 (J2 li) / dnd^dTsdkfdk^dklS (J2 K ) B^(k[,k' 2 , k' 3 ) 



S(k' 1 ,T 1 )S(k' 2 ,T2)S(k' 3 ,T 3 )e- 



ik%(T -Ti) -ik%(T -T 2 ) p -ik%(T -T 3 ) 



where k' is k evaluated with k" = 1/(tq — tr). If we define a new transfer function, 



A (W = J {tq dT _ T? SWm 2 + l 2 Wo - Wrfe*'* , 



(82) 



(83) 



and use, {a(\i)a(h)a(\s)) « (2tt) 2 (5 2 1^) &f^* , then we find that, in the flat sky limit, the expression for the reduced 
bispectrum is, 



rflat Qo - Tr) 

hhh ( 2 7t) 2 



2 /'CXD 



dkfdk z 2 dk!5 (J2 K) B*(k[,k 2 , k' 3 )A(h,kf)A(l 2 , k^)A(l 3 , kl) . 



(84) 
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Figure 11: A comparison of the full calculation of the reduced bispectrum with the flat sky approximation. Top is a comparison 
of bin and we find excellent agreement when / > 150. The bottom two plots show a cross section through h + h + h = 1000 
with the full case on the left and the flat sky approximation on the right. Again we see that the two methods agree except in 
the corners when one of the U < 150. The close agreement between the two independent methods establishes the accuracy of 
the two codes 



For small angles bf 



flat 



bi ± i 2 i 3 and so we can use (84) to calculate the bispectrum. We evaluate the two integrals 



using cubic splines. 

We have calculated the reduced bispectrum using the full method detailed in the previous section, and again in the 
flat sky case, comparing them in figure (11). We find that the flat sky approximation becomes valid when all three 
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h > 150, producing less than 1% error. Also, as the two methods are completely independent of each other, this 
provides a powerful cross check of the accuracy. The flat sky approximation allows us to calculate the bispectrum, 
when all three ^ > 150, more than 300 times faster. This is a dramatic improvement only leaving small subregions 
near the corners and edges which require the full calculation. Nevertheless, even with this much faster method, the 
reduced bispectrum at Planck resolution represents a formidable challenge, unless we significantly reduce the number 
of points at which it needs to be evaluated. 




Figure 12: Transformation to map the tetrahedron to a cube for interpolation. On the left we have the tetrahedron and on the 
right the tetrahedron after the mapping is applied. The blue lines (solid) define the cubic cell over which we interpolate. The 
red point (square) is the value that falls inside the cube after the transformation and is not used for interpolation. Pink points 
and lines (dashed) are the parts of the tetrahedron outside the cell being interpolated. 



Cubic interpolation 

From the plots comparing the full calculation with the flat sky approximation in figure (11), we note that 
G _1 (/i, h^hhh i s actually very smooth. This is expected as for models with smooth shape functions all the 
structure present in the reduced bispectrum must be due to the acoustic peaks in the transfer functions. As a result 
the reduced bispectrum will only contain features that oscillate with periods of I « 200. This indicates that we only 
need calculate the reduced bispectrum on a sparse grid and the remaining points could be generated via interpolation. 
One major problem is in selecting the grid to interpolate over. With the triangle condition on the three U limiting 
us to a tetrahedron we cannot use the usual schemes as when they straddle the boundary they give incorrect results. 
The geometric integral returns zero for I combinations that violate the triangle inequality creating a discontinuity 
which leads to poor convergence. We can circumvent this problem by rotating and stretching the allowed region so it 
forms a rectangular grid via the transform, 

i'i = \{h + h- h) , 

l'i = \ (h + h- h) , 
l'z = \ (h + h- h) ■ 

We can then use cubic interpolation to calculate the remaining points before transforming back to obtain the bis- 
pectrum for all combinations of U. There are some minor issues with this approach. For example, not all points fall 
onto the grid when we rotate. If we are using a grid with steps in I of 10 we would find that the first cell would 
be constructed from: 6222, 3 permutations of 621010, 3 permutations of 6101020, and 6202020- This leaves 6101010 
sitting in the centre of the cell and it is ignored in the subsequent interpolation, see figure (12). As a result we must 
calculate the reduced bispectrum on twice the density of points that we require, but they are not entirely redundant 
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as we use them dynamically as a cross-check of the accuracy of the interpolation. We encounter a similar issue when 
transforming back with the inverse mapping, li = 1' 2 + l f 3 , h = I3 + l[ , h = l[ + Z 2 • 

Interpolation reduces the number of points required to calculate the reduced bispectrum by several orders of 
magnitude and so together with the flat sky approximation and the parallelisation of the code, the problem becomes 
tractable for all the models reviewed in section 3. 




Figure 13: The equal 1 bispectrum for all classes of models. They are, from top to bottom at Z = 220: equilateral, smoothed 
warm, warm, local, smoothed flat, and feature. As all the models scale as k~ 6 they all produce similar bin, with the exception 
of the feature model whose oscillations can be clearly seen. 



CMB BISPECTRUM MODEL COMPARISON 



We have discussed the smooth nature of the CMB bispectrum and its general properties in a previous publication 
[1], so our purpose here is to investigate the distinguishing signatures of the models we outlined in section 3. With the 
method developed in the previous sections we can evaluate the entire reduced bispectrum for any primordial model, 
calculating up to I = 2000 in approximately 100 processor hours; this has been achieved to at least 1% accuracy 
for all the models investigated. In figure 13, we plot the central value bm for all five primordial classes of shapes - 
equilateral (36), local (35), warm (72), flat (73) and feature (74). The first four of these models are scale-invariant, 
so the bm all take broadly the same profile but with different normalisations. We note that this figure demonstrates 
the oscillatory properties of the transfer functions which, as for the CMB power spectrum, create a series of acoustic 
peaks around I = 200, 500, 800, .... There is a stark contrast, however, with the feature model which has a non-trvial 
scaling. Initially, it is anticorrelated with the other shapes, so that the primary peak has opposite sign, however, for 
increasing Z the phase of the oscillations becomes positively correlated by the second and third peaks (this, of course, 
reflects the particular choice of k* in (74)). 

Of course, to observe the key differences between the models we must study the bispectrum in the plane orthogonal 
to the (Z, Z, indirection, that is, the directions reflecting changes in the primordial shape functions. To this end, 
in figures 14 and 15 we plot cross-sectional slices through the reduced bispectrum; we choose triangles satisfying 
^1+^2 + ^3 = 1000 just beyond the primary peak so as to reduce the effect of the transfer functions relative to the 
primordial shape. For these slices and in the subsequent figures 16-18 for the full three-dimensional bispectrum, we 
divide bm by D(li, Z2, Z3), the large-angle CMB bispectrum solution (26) for a primordial shape which is constant. This 
is analogous to multiplying the power spectrum C/'s by 1(1 + 1), because it serves to flatten the bispectrum except for 
the effect of the non-constant primordial shape (and the oscillating transfer functions). (We note that in our previous 
paper [1], for plotting purposes we divided all the CMB bispectra by the large-angle local solution G(h, Z2, h) given 
in (20), but this is not as useful for comparison purposes.) 

The bispectrum slices for the different models shown in figure (14) directly mimic the primordial shapes from which 
they originated. The equilateral model has the majority of its signal under a prominent central peak (i.e. equilateral 
triangles), whereas the local model has a nearly flat interior with most signal at sharp peaks in the corners (i.e. for 
squeezed triangles). As expected, the (smoothed) flat model is strongly peaked along the edges for flattened triangle 
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configurations. The last model illustrated in figure (14) is the feature model whose primordial shape function is 
constant across k = const, slices. On I = const, slices, therefore, the feature model should behave like the constant 
model, showing only the effect of the transfer functions in its structure. Finally, there is a slice through the warm 
model shown in figure (15) in both original and smoothed versions. These have signal when the ^ are in several 
configurations, squeezed, flattened and equilateral. Note the strong effect that smoothing has on suppressing the 
dominant contribution from squeezed states. A better understanding of the asymptotic behaviour in squeezed and 
flattened limits is clearly necessary if we are to make robust quantitative predictions for both the warm and fiat 
models. 

To understand the differences between the four models in the same equilateral class, it is easiest to look at the 3D 
plots of the reduced bispectrum shown in figure 16. For all models we see the same peak structure with the maximums 
at positions where all the three U correspond to peaks in the power spectrum. The largest peak is then when all three 
li = 220, corresponding to the large blue region near the origin. The four models are most strongly differentiated 
by their behaviour in the flattened limit, as can be seen in figure (16). As a result they can be separated by their 
behaviour around the peak positioned where l\ = I2 = 250, I3 = 500, (the magenta part just above the U = 220 
peak). In the equilateral case it is a modest feature but in DBI and single it is larger, connecting up to create a forked 
range of structures near the faces of the tetrahedron. For ghost inflation, which becomes negative as we approach the 
flattened limit, the peak is almost non-existent. 

Despite the visual differences between equilateral models, we note that the CMB correlator (42) accurately confirms 
the strong correlations forecast by our simple shape correlator (52), see table (II). As predicted, all three shapes - 
DBI, ghost and single - correlate closely with the phenomenological equilateral shape to within 96% (this agrees with 





Figure 14: The reduced bispectrum for different classes of models plotted on slices where h + h + ^3 = 1000. Clockwise from 
top left, Local, Equilateral, FlatS, and Feature. The effect of the primordial shape function can be clearly seen in the resulting 
bispectra. 
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Figure 15: The reduced bispectrum for warm (left) and warmS (right) models plotted on slices where l\ + I2 + ^3 = 1000. 

[13] for the DBI model, the only direct calculation elsewhere of a non-separable case). Again the largest difference 
is between the ghost and single models with an 89% correlation, which was slightly under-estimated by the shape 
correlator at 85%, presumably because of subtle changes in weighting arising from the transfer functions, as well 
as smoothing effects. We confirm that these four equilateral shapes will be very difficult to differentiate without a 
bispectrum detection of high significance. 

The local model is shown in figure (17), demonstrating a marked constrast with the 3D equilateral bispectrum. 
The dominance of signal in the squeezed limit creates strong parallel ridges which connect up and emanate along the 
corner edges of the tetrahedron. The 51% CMB correlation between the local and equilateral models is reasonably 
well predicted by the shape correlator at 41%, though again underestimated. The full results for all the CMB cross- 
correlators can be seen in Table II. (We note that the estimated 2D cosines between the local model and the DBI 
and ghost models are in qualitative agreement with ref. [11], see also [13].) 

It remains to briefly discuss the bispectrum and correlations of the last three models. The two 3D plots for warm 
inflation, see figure (18), seem superficially similar but the effect of smoothing is to move signal from the squeezed 
states into the centre. This apparently minor shift is enough to reduce the fisher correlation between the warm and 
warmS to 48%, providing the most significant outlier for the shape correlator with an 80% correlation. It demonstrates 
the sensitivity of the analysis to the exact cut-off used for the primordial bispectrum in the squeezed limit, but it may 
also be accentuated by an increased weight around the edges implicit in the CMB estimator (see earlier discussion in 
section 3. The smoothed flat model bispectrum shown in (17) has become visually more similar to the local bispectrum 
than might be expected from the primordial shapes, with the CMB correlator reflecting this change at 79%. The flat 
and local shapes could be difficult to disentangle, though we note that this result depends on the imposition of an 
arbitrary cut-off to regularise the flat shape. 

Finally, the feature model clearly represents an entirely distinct type of bispectrum, which is evident from its very 
different behavior in the (1,1, /)-direction. The anticorrelation of the primary peak, relative to the other peaks, is clear 
in figure (19) from the nodal plane which cuts across the bispectrum. The poor correlation predicted with all the 
other shapes, is confirmed by the CMB cross-correlations which are all below 45%. (Here, the interplay with the nodal 
points introduced by the transfer functions makes these results strongly dependent on l m ax)- Preliminary analysis 
based on the approximate analytic shapes for an oscillatory model indicates a further independent shape which could 
be distinguished from the other classes given a reasonably significant bispectrum detection. 

In figure 20, we have plotted the full CMB bispectrum correlator (42) against the simple shape correlator (52) for all 
the models investigated, demonstrating their remarkable concordance; highly correlated shapes agree accurately, while 
the shape correlator understimates the correlation of independent shapes (usually by about 5-15%). Such a simple 
predictor of model correlations is important given the computational effort required to compare the CMB bispectra 
directly. This analysis confirms that there are indeed five distinct classes of models among the cases reviewed: The 
equilateral class, the local class (which includes the non-local model), the warm model, the flat class and the feature 
model. 

Here we note the importance of the weight, eqn. (51), in the shape correlator. For the feature model it is vital to 
achieve the correct correlation with any other choice producing extremely poor results. For the scale invariant models 
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Figure 16: 3D Plots of the reduced bispectrum for models in the equilateral class. Here, and in the following, we have plotted 
surfaces of equal bispectra after division by the constant analytic solution. Clockwise from top left we have equilateral, DBI, 
Ghost and Single. 

it is not as important to have the correct weight, as long as it weights each point on a slice equally. However, while 
tolerable results can still be achieved for the scale invariant models with weights uo = 1 and uj = l/(&i + + fe) 2 , we 
still see much better agreement between the two correlators when we use the correct weight, uo = l/(&i + &2 + fe), 
due to the effect of the shape of the region of integration as discussed earlier. 



27 




Figure 17: 3D Plots of the reduced bispectrum for the local model (left) and the flattened model after smoothing (right). 




Figure 18: 3D Plots of the reduced bispectrum for warm inflation before smoothing, on the left, and after smoothing, on the 
right. 

CMB BISPECTRUM NORMALIZATION 

The question now is how best to compare and contrast observational limits for such a wide variety of possible 
models. We need a framework for deciding a sensible definition of /nl for each class of models, normalised so that 
comparisons can be made. There have been many attempts at extending /jy£ ai to other models in the literature and 
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Table II: Comaprison of bispectrum and shape correlators 



Modell- Model2 


C(S,S') 


C{B,B') 


Modell- 


Model2 


C(S,S') 


C(B,B') 


DBI - Equi 


0.99 


0.99 


Warm - 


DBI 


0.38 


0.39 


DBI - Flat S 


0.39 


0.48 


Warm - 


Equi 


0.44 


0.42 


DBI - Ghost 


0.94 


0.95 


Warm - 


Flat S 


0.01 


0.21 


DBI - Local 


0.50 


0.56 


Warm - 


Ghost 


0.50 


0.43 


DBI - Single 


0.98 


0.99 


Warm - 


Local 


0.30 


0.52 


DBI - Warm S 


0.55 


0.69 


Warm - 


Single 


0.29 


0.35 


Equi - Flat S 


0.30 


0.39 


Warm - 


Warm S 


0.80 


0.48 


Equi - Ghost 


0.98 


0.98 










Equi - Local 


0.46 


0.51 


Feature 


- DBI 


-0.41 


-0.44 


Equi - Single 


0.95 


0.96 


Feature 


- Equi 


-0.36 


-0.43 


Equi - Warm S 


0.63 


0.76 


Feature 


- Flat S 


-0.44 


-0.32 


Flat S - Ghost 


0.15 


0.24 


Feature 


- Ghost 


-0.26 


-0.42 


Flat S - Local 


0.62 


0.79 


Feature 


- Local 


-0.41 


-0.39 


Flat S - Single 


0.49 


0.60 


Feature 


- Single 


-0.46 


-0.44 


Flat S - Warm S 


-0.03 


-0.04 


Feature 


- Warm 


-0.05 


-0.27 


Ghost - Local 


0.37 


0.42 


Feature 


- Warm S 


-0.08 


-0.20 


Ghost - Single 


0.86 


0.89 










Ghost - Warm S 


0.71 


0.82 










Local - Single 


0.55 


0.62 










Local - Warm S 


0.27 


0.14 










Single - Warm S 


0.44 


0.58 











the one most extensively used is to normalise the shape functions against a single point. This is done for the equilateral 
model and for the warm model, (although for the warm model, as they work with the curvature perturbation £, rather 
than <£, there is an additional factor of 3/5 in their definition). Creating a new non-Gaussianity parameter for each 
model, fx™ and /^J™, the models are normalised such that, 

S local (k, k, k) = S equi (k, k, k) = ls warm (k, k, k) . (85) 

o 

As all three shapes have the same scaling, this ratio is independent of fc, and the method has the benefit of simplicity. 
However, in the equilateral model we are normalising the maximum of one shape to the minimum of the other (local), 
so we find that a similar level of non-Gaussianity produces an f^f™ that is over three times larger. For the warm 
case, we have a /$£ rm that is almost four times larger, than /^£ aZ . Consequently bounds on and /^^ rm seem 

far weaker than those for /jy£ ai which is entirely misleading, see eqns (1-3). The approach makes even less sense when 
the shapes have a running, in which case an arbitrary value for k must be used for the normalisation, and there is no 
obvious extension to models with features or oscillatory behaviour. 

One approach to normalisation would be to define Jnl such that a given model with /atl = 1 produces the same 
level of bispectrum signal as the local case, also with f 1 ^ 1 = 1, i.e. 

M{B) =N{B local ), (86) 

where M is denned by equation (28) which we repeat here for convenience, 



\ y ^hhh (87) 



We could regard M as the total integrated bispectrum signal of a particular model. Note that here we calculate J\f for 
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Figure 19: 3D Plot of the reduced bispectrum for the feature model. 



an ideal experiment without noise or beam effects covering the full I range of the non-Gaussian signal; this would then 
give error bars which are approximately the same for all models. As it is defined at late times, this would also have 
the advantage that it could be used for the bispectrum from cosmic strings, or other general phenomena, e.g. allowing 
an effective f^™ 9 to be defined and related to G/i, the string tension. 

When we calculate the estimator from equation (27) for a particular experiment, including its particular noise 
and beam profiles, we will find any loss of sensitivity to the bispectrum for a particular model being reflected by an 
increase in the error bars. Thus the errors on models are now directly related to the sensitivity of the experiment to 
that model, rather than to an arbitrary choice for the definition of f^if el . 

To avoid calculating the full bispectrum for each model it remains more straight forward to define the normalisation 
relative to the level of primordial signal instead. Given the success of the correlator based on the shape function (52), 
we can use a normalisation based on the integrated non-Gaussianity in /c-space. As S(ki : k2 : ks) represents the 
primordial signal that we project forward to obtain the CMB bispectrum today, then using the same scaling as in the 
shape correlator should represent a more consistent measure of the actual non-Gaussianity measured in the CMB, 
that is, 



M 2 (S) = / S 2 (k u k 2 ,k 3 )w(k u k 2 ,k 3 )dV k . 



(88) 



where Vk is the allowed region in fc-space with the fcj < k max (see figure 1) and we take the weight function 
w(k\,k2,ks) = (k\ + ki + k^) -1 , as before replicating the ^-dependence of the CMB estimator (50). (A similar 
suggestion, but with a different weight, was made in ref. [12].) We then propose to normalise the /nls in each class 
of models relative to the local model, that is, such that 



Af(S) = Af(S iocai ). 



We find then that the 2<r limits quoted above become, 



-4 < m < 80 , 

-34</X<57, 
-107 < f% a L rm < 11. 



(89) 

(90) 
(91) 
(92) 
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Figure 20: Comparison of the bispectrum correlator and primordial shape correlator. The blue line represents perfect agreement 
between the two, with the red data points (circles) showing correlations between: DBI, Equi, FlatS, Ghost, Local, Single, and 
WarmS. The correlators are in good agreement in these cases. The green data points (triangles) are correlations between these 

7 models and the unsmoothed warm shape, while the purple data points (squares) are the correlations between the previous 

8 models and the feature model. For the unsmoothed warm model, there is one significant outlier when the pathological 
behaviour in the squeezed limit for the warm case is being picked up in the shape correlator but not in the bispectrum (it has 
been smoothed by projection). Warm and WarmS models are poorly correlated indicating that the squeezed limit must be 
understood fully before we can consider any constraints to be robust. The feature model is anti-correlated will all other models 
at about the 40% level or below. 



where the standard deviation is now 21, 23, and 29 respectively and we use Jnl to denote that is it denned using the 
new approach to normalisation. We can see that the local and equilateral models are constrained at the same level 
with the warm constraint being weaker as the estimator used in [3] is not optimal in the presence of inhomogeneous 
noise. This normalisation is applicable to all models regardless of their form, it is simple to calculate, and the Jnl 
for each class of models then represents a similar level of the measurable non-Gaussian signal. 

Equally, we can use this standardised normalisation, together with the correlator results in table II, to naively 
forecast Jnl and its errors for alternative models which are not yet constrained. Supposing our universe actually 
possessed significant local non-Gaussianity, then on the basis of local estimator observations with Jnl — 38 ±21 as in 
(90), the equilateral model should yield ~ 17 ±21 (consistent with the observed result), while for the flat model 

« 24 ± 21. Conversely, if our universe possessed flat non-Gaussianity, then given the local result (90) we might 
obtain a marginally significant result ff lat ^ 61 ± 21. We conclude that discovery potential remains, even with the 
present CMB data, for the independent shapes which have not yet been fully investigated, such as feature models 1 . 



We note that the warm inflation results do not seem to match expectations from the shape correlator, suggesting an anticorrelation 
with the local model, rather than the 30% correlation of table II. However, we caution that the warm result has not been independently 
verified and it also depends sensitively on arbitrary cut-offs imposed on the shape function. We note that we have removed the factor 
of 3/5 from the definition of f%l rm in 90. 
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To reiterate the value of a standardised normalisation (88), we note an obvious failing of the previous method. This 
comes from extending an observational result for one model to all the highly correlated models in the same class, 
such as equilateral to DBI, ghost and single. If we normalise the models using the centre point S(fc, fc), as is done 
currently in the literature, then with the different N(S) in each case the original equilateral limit f^™ = 51 ± 101 
transfers to the following inconsistent limits 

/££ 7 = 47±93, (93) 
/£2' = 40±78, (94) 
f g N h ° st = 59 ±116. (95) 

In contrast, if we normalise using N(S) in (88), then the equilateral limit (91) on Jnl transfers across nearly identically 
in all these cases, because of their 95% cross-correlation. 



CONCLUSIONS 



In this paper we have presented a comprehensive set of formalisms for comparing, evolving, and constraining 
primordial non-Gaussian models through the CMB bispectrum. 

The primary goal was to directly calculate CMB bispectra from a general primordial model, enhancing methods 
previously outlined in ref. [1]. This was achieved by exploiting common features of primordial bispectra to reduce the 
dimensionality of the transfer functions required to evaluate the CMB bispectra. The new innovations reported here 
include the use of the flat sky approximation when all l{ > 150, greatly reducing computational times for most of 
the allowed region, and a cubic reparametrisation, significantly reducing the number of points required for accurate 
interpolation of the bispectrum. (We note that this CMB bispectrum code is being prepared for public distribution.) 
These methods make feasible the repetitive calculation of highly accurate CMB bispectra at Planck resolution without 
specific assumptions about the separability of the underlying primordial bispectrum. 

Further, we have calculated the CMB bispectra for all the distinct primordial bispectrum shapes S(ki, &2, ks) 
currently presented in the literature, motivated by a range of inflationary and other cosmological scenarios. We have 
presented these, plotted relative to the large- angle CMB bispectrum for the constant primordial model (S const - = 1) 
for which we presented an analytic solution (25). The CMB bispectra from the different primordial models exhibit a 
close correspondence to their original shape modulated, of course, by the oscillatory transfer functions. 

We were able to quantitatively determine the observational independence of the CMB bispectra by measuring their 
cross-correlations using the estimator (27). These results revealed five independent classes of shapes which it should 
be possible to distinguish from one another with a significant detection of non-Gaussianity in future experiments such 
as Planck. These were the equilateral (36), local (35), warm (72), and flat (73) shapes, all described in section 3, 
together with the feature model (74) which is basically the constant shape (21) plus broken scale-invariance. Different 
models belonging within the same class will be difficult to distinguish, a fact best demonstrated by the 95% correlation 
of the equilateral shape with DBI (68), ghost (69) and single (71) shapes. 

We have also presented a shape correlator (52) which provides a fast and simple method for determining the 
independence of different shapes. In particular, for highly correlated models, it yielded results accurately reflecting 
those of the full CMB correlator (42), thus avoiding unnecessary calculation. The shape correlator also reliably 
identified poorly correlated models, that is, new shapes for which a full CMB bispectrum analysis was warranted. 
We also proposed a straightforward two-dimensional eigenmode analysis of shape functions, valid for nearly scale- 
invariant models. This allowed us to identify the shape correlations with products of eigenvalues of the non-orthogonal 
eigenmodes, immediately revealing, for example, why warm and local models are independent. In principle, the 
analysis can guide the theoretical search for primordial models with distinctive non-Gaussian signatures. It also 
revealed the constant eigenmode (or shape (21)) as the primary cause of the cross-correlation between many models, 
suggesting strategies for distinguishing, for example, between local and equilateral models. 

Finally, we proposed an alternative approach (88) to the normalisation of the non-Gaussianity parameter f^L using 
the shape autocorrelator. This new normalisation allows us to employ Jnl to systematically compare the true level 
of non-Gaussianity in different models. In contrast, with current methods, the constraints for competing models 
can vary by a factor of four or more, with bounds varying significantly for models even in the same class of highly 
correlated shapes. 

A detection of non-Gaussianity would have profound consequences for our understanding of the early universe, 
uprooting the present simplest inflationary paradigm. The present work indicates that the next generation of CMB 
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experiments (notably Planck) may be able to distinguish between different classes of shapes for primordial non- 
Gaussianity. Delineating the bispectrum shape would provide important clues about viable alternative theories for 
the origin of large-scale structure in the universe. 
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